! ! This writes output from simulations
module simwrite

use params
use simulate

implicit none

integer, private :: i, j, k, t, c1, c2, c3, c4, c5, k1, k2, k3, k4, k5
real(prec), private :: nalive,nalivec(ncohort),nalivectot(ncohort),nFTcohort(ncohort),nmarried,allunhealthy,allalive

contains

!------------------
subroutine writesim

integer outcsv
integer nsimc(ncohort)
double precision   :: outmat(nt+1,73)
double precision, allocatable :: alivewlth(:,:),alivewlthc1(:,:),alivewlthc2(:,:),alivewlthc3(:,:),alivewlthc4(:,:),alivewlthc5(:,:)
double precision, allocatable :: assk_c1(:),assk_c2(:),assk_c3(:),assk_c4(:),assk_c5(:)
double precision, allocatable :: alivemarr(:),aliveconp(:),aliveconk(:)
double precision :: practice(10), prac25, prac50,prac75
double precision :: pct10,pct25,pct50,pct75,pct90

! averages, medians from simulation output
avinfoverall=0
avFToverall=0
avFTsickoverall=0
avFTwelloverall=0
avPToverall=0
avPTsickoverall=0
avPTwelloverall=0
avSIpoverall=0
avwealthp50=0
avwealthk50=0
allalive=0
allunhealthy=0
nalivectot=0
nFTcohort=0
do t = 1,nt
    nalive = dble(totalive(t))    !divide by number living
	do i=1,ncohort
		nalivec(i)= dble(totalivec(t,i))
	enddo	
	if (t>=1 .and. t<=7) then
		nalivectot(1)= nalivectot(1)+dble(totalivec(t,1))
		nFTcohort(1)= nFTcohort(1)+dble(totFTcohort(t,1))
	endif
	if (t>=4 .and. t<=10) then
		nalivectot(2)= nalivectot(2)+dble(totalivec(t,2))
		nFTcohort(2)= nFTcohort(2)+dble(totFTcohort(t,2))
	endif
	if (t>=6 .and. t<=12) then 
		nalivectot(3)= nalivectot(3)+dble(totalivec(t,3))
		nFTcohort(3)= nFTcohort(3)+dble(totFTcohort(t,3))
	endif	
	if (t>=9 .and. t<=15) then
		nalivectot(4)= nalivectot(4)+dble(totalivec(t,4))
		nFTcohort(4)= nFTcohort(4)+dble(totFTcohort(t,4))
	endif	
	if (t>=11) then
		nalivectot(5)= nalivectot(5)+dble(totalivec(t,5))
		nFTcohort(5)= nFTcohort(5)+dble(totFTcohort(t,5))
	endif		
	
	
	nmarried = dble(totmar(t))    !divide by number living
    ! Decision variables
    avvalp(t) = totvalp(t)/nalive
    avvalk(t) = totvalk(t)/nalive
    avconp(t) = totconp(t)/nalive
    avconk(t) = totconk(t)/nalive
    avincp(t) = totincp(t)/nalive
    avinck(t) = totinck(t)/nosims    
    avsavp(t) = totsavp(t)/nalive
    avsavk(t) = totsavk(t)/nosims
    avSIp(t)  = totSIp(t)/nalive
    avSIk(t)  = totSIk(t)/nosims
    avhlthy(t)= tothlthy(t)/nalive
    avmar(t)  = totmar(t)/nalive
    avmu(t)   = totmu(t)/nmarried
    avkap(t)  = totkap(t)/nmarried
    avPT(t)   = totPT(t)/nalive
    avFT(t)   = totFT(t)/nalive

    avFT_sick(t) = totFT_sick(t)/(dble(totalive(t))-tothlthy(t))
    avFT_well(t) = totFT_well(t)/dble(tothlthy(t))
	
    avPT_sick(t) = totPT_sick(t)/(dble(totalive(t))-tothlthy(t))
    avPT_well(t) = totPT_well(t)/dble(tothlthy(t))	
	
    !median -- must take out dead people for parents! also do by cohort! but don't take out dead people for kids...
    allocate(alivewlth(int(nalive),2))
    alivewlth=0
	allocate(alivewlthc1(int(nalivec(1)),2))
	alivewlthc1=0
	allocate(alivewlthc2(int(nalivec(2)),2))
	alivewlthc2=0
	allocate(alivewlthc3(int(nalivec(3)),2))
	alivewlthc3=0
	allocate(alivewlthc4(int(nalivec(4)),2))
	alivewlthc4=0
	allocate(alivewlthc5(int(nalivec(5)),2))
	alivewlthc5=0
	allocate(assk_c1(int(totcohort(1))))
	assk_c1=0
	allocate(assk_c2(int(totcohort(2))))
	assk_c2=0
	allocate(assk_c3(int(totcohort(3))))
	assk_c3=0
	allocate(assk_c4(int(totcohort(4))))
	assk_c4=0
	allocate(assk_c5(int(totcohort(5))))
	assk_c5=0	
    allocate(aliveconp(int(nalive)))
    aliveconp=0
    allocate(aliveconk(int(nalive)))
    aliveconk=0
    allocate(alivemarr(int(nmarried)))
    alivemarr=0
    j=1
    k=1
	c1=1
	c2=1
	c3=1
	c4=1
	c5=1
	k1=1
	k2=1
	k3=1
	k4=1
	k5=1
	nsimc=0
    do i=1,nosims
        if (alivemat(i,t)==1) then
            alivewlth(j,1)=assp(i,t) 
            alivewlth(j,2)=assk(i,t) 
            aliveconp(j)=conp(i,t)
            aliveconk(j)=conk(i,t)
            j=j+1
			if (simcohort(i)==1) then
				alivewlthc1(c1,1)=assp(i,t)
				alivewlthc1(c1,2)=assk(i,t)
				c1=c1+1
			elseif (simcohort(i)==2) then
				alivewlthc2(c2,1)=assp(i,t)
				alivewlthc2(c2,2)=assk(i,t)
				c2=c2+1	
			elseif (simcohort(i)==3) then
				alivewlthc3(c3,1)=assp(i,t)
				alivewlthc3(c3,2)=assk(i,t)
				c3=c3+1					
			elseif (simcohort(i)==4) then
				alivewlthc4(c4,1)=assp(i,t)
				alivewlthc4(c4,2)=assk(i,t)
				c4=c4+1		
			elseif (simcohort(i)==5) then
				alivewlthc5(c5,1)=assp(i,t)
				alivewlthc5(c5,2)=assk(i,t)
				c5=c5+1	
			endif
        endif
		if (simcohort(i)==1) then
			nsimc(1)=nsimc(1)+1
			assk_c1(k1)=assk(i,t)
			k1=k1+1
		elseif (simcohort(i)==2) then
			nsimc(2)=nsimc(2)+1
			assk_c2(k2)=assk(i,t)
			k2=k2+1
		elseif (simcohort(i)==3) then
			nsimc(3)=nsimc(3)+1
			assk_c3(k3)=assk(i,t)
			k3=k3+1	
		elseif (simcohort(i)==4) then
			nsimc(4)=nsimc(4)+1
			assk_c4(k4)=assk(i,t)
			k4=k4+1
		elseif (simcohort(i)==5) then
			nsimc(5)=nsimc(5)+1
			assk_c5(k5)=assk(i,t)
			k5=k5+1			
		endif
		if (marstat(i,t)==1) then
            alivemarr(k)=kappct(i,t)
            k=k+1
        endif
    enddo
    call PERCENTILES(alivewlth(:,1),int(nalive),pct10,pct25,pct50,pct75,pct90)
    wealthp10(t)=pct10
    wealthp25(t)=pct25
    wealthp50(t)=pct50
    wealthp75(t)=pct75
    wealthp90(t)=pct90
	if (t>=1 .and. t<=9) then
		call PERCENTILES(alivewlthc1(:,1),int(nalivec(1)),pct10,pct25,pct50,pct75,pct90)
		wealthcp10(1,t)=pct10
		wealthcp25(1,t)=pct25
		wealthcp50(1,t)=pct50
		wealthcp75(1,t)=pct75
		wealthcp90(1,t)=pct90
	endif
	if (t>=4 .and. t<=12) then
		call PERCENTILES(alivewlthc2(:,1),int(nalivec(2)),pct10,pct25,pct50,pct75,pct90)
		wealthcp10(2,t-3)=pct10
		wealthcp25(2,t-3)=pct25
		wealthcp50(2,t-3)=pct50
		wealthcp75(2,t-3)=pct75
		wealthcp90(2,t-3)=pct90
	endif
	if (t>=6 .and. t<=14) then 
		call PERCENTILES(alivewlthc3(:,1),int(nalivec(3)),pct10,pct25,pct50,pct75,pct90)
		wealthcp10(3,t-5)=pct10
		wealthcp25(3,t-5)=pct25
		wealthcp50(3,t-5)=pct50
		wealthcp75(3,t-5)=pct75
		wealthcp90(3,t-5)=pct90
	endif
	if (t>=9 .and. t<=15) then
		call PERCENTILES(alivewlthc4(:,1),int(nalivec(4)),pct10,pct25,pct50,pct75,pct90)
		wealthcp10(4,t-8)=pct10
		wealthcp25(4,t-8)=pct25
		wealthcp50(4,t-8)=pct50
		wealthcp75(4,t-8)=pct75
		wealthcp90(4,t-8)=pct90
	endif
	if (t>=11 .and. t<=18) then
		call PERCENTILES(alivewlthc5(:,1),int(nalivec(5)),pct10,pct25,pct50,pct75,pct90)
		wealthcp10(5,t-10)=pct10
		wealthcp25(5,t-10)=pct25
		wealthcp50(5,t-10)=pct50
		wealthcp75(5,t-10)=pct75
		wealthcp90(5,t-10)=pct90
	endif
	
   	if (t>=1 .and. t<=9) then
		call PERCENTILES(assk_c1(:),totcohort(1),pct10,pct25,pct50,pct75,pct90)
		wealthck10(1,t)=pct10
		wealthck25(1,t)=pct25
		wealthck50(1,t)=pct50
		wealthck75(1,t)=pct75
		wealthck90(1,t)=pct90
	endif
	if (t>=4 .and. t<=12) then
		call PERCENTILES(assk_c2(:),totcohort(2),pct10,pct25,pct50,pct75,pct90)
		wealthck10(2,t-3)=pct10
		wealthck25(2,t-3)=pct25
		wealthck50(2,t-3)=pct50
		wealthck75(2,t-3)=pct75
		wealthck90(2,t-3)=pct90
	endif
	if (t>=6 .and. t<=14) then 
		call PERCENTILES(assk_c3(:),totcohort(3),pct10,pct25,pct50,pct75,pct90)
		wealthck10(3,t-5)=pct10
		wealthck25(3,t-5)=pct25
		wealthck50(3,t-5)=pct50
		wealthck75(3,t-5)=pct75
		wealthck90(3,t-5)=pct90
	endif
	if (t>=9 .and. t<=15) then
		call PERCENTILES(assk_c4(:),totcohort(4),pct10,pct25,pct50,pct75,pct90)
		wealthck10(4,t-8)=pct10
		wealthck25(4,t-8)=pct25
		wealthck50(4,t-8)=pct50
		wealthck75(4,t-8)=pct75
		wealthck90(4,t-8)=pct90
	endif
	if (t>=11 .and. t<=18) then
		call PERCENTILES(assk_c5(:),totcohort(5),pct10,pct25,pct50,pct75,pct90)
		wealthck10(5,t-10)=pct10
		wealthck25(5,t-10)=pct25
		wealthck50(5,t-10)=pct50
		wealthck75(5,t-10)=pct75
		wealthck90(5,t-10)=pct90
	endif	
	
    call PERCENTILES(aliveconp(:),int(nalive),pct10,pct25,pct50,pct75,pct90)
    medconp(t)=pct50  
    call PERCENTILES(aliveconk(:),int(nalive),pct10,pct25,pct50,pct75,pct90)
    medconk(t)=pct50  
    deallocate(alivewlth)
	deallocate(alivewlthc1)
	deallocate(alivewlthc2)
	deallocate(alivewlthc3)
	deallocate(alivewlthc4)
	deallocate(alivewlthc5)
	deallocate(assk_c1)
	deallocate(assk_c2)
	deallocate(assk_c3)
	deallocate(assk_c4)
	deallocate(assk_c5)	
	deallocate(aliveconp)
    deallocate(aliveconk)
    if (nmarried>=1) then 
        call PERCENTILES(alivemarr(:),int(nmarried),pct10,pct25,pct50,pct75,pct90)
        kappap25(t)=pct25
        kappap50(t)=pct50
        kappap75(t)=pct75
    endif

    do k=1,5
        avins(k) = totins(1,k)/totins(2,k) 
		avinsPO(k)= totinsPO(1,k)/totinsPO(2,k) 
    enddo
    avins(6)= (totins(1,1)+totins(1,2)+totins(1,3)+totins(1,4)+totins(1,5))/dble(totalive(1))
    avinsPO(6)= (totinsPO(1,1)+totinsPO(1,2)+totinsPO(1,3)+totinsPO(1,4)+totinsPO(1,5))/dble(totalive(1))

    do k=1,5
        avinf(k) = totinf(1,k)/totinf(2,k)
    enddo
    avinf(6)=(totinf(1,1)+totinf(1,2)+totinf(1,3)+totinf(1,4)+totinf(1,5)) / &
             (totinf(2,1)+totinf(2,2)+totinf(2,3)+totinf(2,4)+totinf(2,5))

    deallocate(alivemarr)

    avFToverall=avFToverall+totFT(t)
	avFTsickoverall=avFTsickoverall+totFT_sick(t)
	avFTwelloverall=avFTwelloverall+totFT_well(t)
	avPToverall=avPToverall+totPT(t)
	avPTsickoverall=avPTsickoverall+totPT_sick(t)
	avPTwelloverall=avPTwelloverall+totPT_well(t)	
    avSIpoverall=avSIpoverall+totSIp(t)

    allalive=allalive+nalive
    allunhealthy=allunhealthy+(nalive-tothlthy(t))

    avwealthp50=avwealthp50+wealthp50(t)
    avwealthk50=avwealthk50+wealthk50(t)


enddo       ! End loop over time

avinfoverall=avinfoverall/allunhealthy
avFToverall=avFToverall/allalive
avFT(nt+1)=avFToverall
avFTsickoverall=avFTsickoverall/allunhealthy
avFT_sick(nt+1)=avFTsickoverall
avFTwelloverall=avFTwelloverall/(allalive-allunhealthy)
avFT_well(nt+1)=avFTwelloverall
avPToverall=avPToverall/allalive
avPT(nt+1)=avPToverall
avPTsickoverall=avPTsickoverall/allunhealthy
avPT_sick(nt+1)=avPTsickoverall
avPTwelloverall=avPTwelloverall/(allalive-allunhealthy)
avPT_well(nt+1)=avPTwelloverall
avSIpoverall=avSIpoverall/allalive
avSIp(nt+1)=avSIpoverall
avwealthp50=avwealthp50/real(nt)
avwealthk50=avwealthk50/real(nt)

do i=1,ncohort
	avFTcohort(i) = nFTcohort(i)/nalivectot(i)
end do

if (pid==0 .and. DoEstimation==0) then 
    !Create output matrix
    do t=1,nt
        outmat(t,1)=t           !age
    end do
    outmat(1:nt,2)=medconp(:)       !mean consumption by age for parents
    outmat(1:nt,3)=medconk(:)       !mean consumption by age for kids
    outmat(1:nt,4)=avincp(:)        !mean income by age for parents
    outmat(1:nt,5)=avinck(:)        !mean income by age for kids
    outmat(1:nt+1,6)=avFT(:)       !mean labor supply
    outmat(1:nt+1,7)=avFT_sick(:)  !mean labor supply
    outmat(1:nt+1,8)=avFT_well(:)  !mean labor supply
	outmat(1:nt+1,9)=avPT(:)
	outmat(1:nt+1,10)=avPT_sick(:)
	outmat(1:nt+1,11)=avPT_well(:)
    outmat(1:6,12)=avinf(:)         !mean informal care by initial asset quintile (unhealthy only)
    outmat(1:6,13)=avins(:)         !mean insurance for low income parents
    outmat(1:nt,14)=avmar(:)     	!mean marital status
    outmat(1:nt,15)=avmu(:)   		!mean pareto weight
    outmat(1:nt+1,16)=avSIp(:)      !mean medicaid parent
    outmat(1:nt+1,17)=avSIk(:)      !mean medicaid kid
    outmat(1:nt,18)=avhlthy(:)
    outmat(1:nt,19)=kappap25(:)
    outmat(1:nt,20)=kappap50(:)
    outmat(1:nt,21)=kappap75(:)
	outmat(1:9,22)=wealthcp10(1,:)
	outmat(1:9,23)=wealthcp25(1,:)
	outmat(1:9,24)=wealthcp50(1,:)
	outmat(1:9,25)=wealthcp75(1,:)
	outmat(1:9,26)=wealthcp90(1,:)
	outmat(1:9,27)=wealthcp10(2,:)
	outmat(1:9,28)=wealthcp25(2,:)
	outmat(1:9,29)=wealthcp50(2,:)
	outmat(1:9,30)=wealthcp75(2,:)	
	outmat(1:9,31)=wealthcp90(2,:)	
	outmat(1:9,32)=wealthcp10(3,:)
	outmat(1:9,33)=wealthcp25(3,:)
	outmat(1:9,34)=wealthcp50(3,:)
	outmat(1:9,35)=wealthcp75(3,:)	
	outmat(1:9,36)=wealthcp90(3,:)	
	outmat(1:9,37)=wealthcp10(4,:)
	outmat(1:9,38)=wealthcp25(4,:)
	outmat(1:9,39)=wealthcp50(4,:)
	outmat(1:9,40)=wealthcp75(4,:)	
	outmat(1:9,41)=wealthcp90(4,:)	
	outmat(1:9,42)=wealthcp10(5,:)
	outmat(1:9,43)=wealthcp25(5,:)
	outmat(1:9,44)=wealthcp50(5,:)
	outmat(1:9,45)=wealthcp75(5,:)	
	outmat(1:9,46)=wealthcp90(5,:)	
	outmat(1:9,47)=wealthck10(1,:)
	outmat(1:9,48)=wealthck25(1,:)
	outmat(1:9,49)=wealthck50(1,:)
	outmat(1:9,50)=wealthck75(1,:)
	outmat(1:9,51)=wealthck90(1,:)
	outmat(1:9,52)=wealthck10(2,:)
	outmat(1:9,53)=wealthck25(2,:)
	outmat(1:9,54)=wealthck50(2,:)
	outmat(1:9,55)=wealthck75(2,:)	
	outmat(1:9,56)=wealthck90(2,:)	
	outmat(1:9,57)=wealthck10(3,:)
	outmat(1:9,58)=wealthck25(3,:)
	outmat(1:9,59)=wealthck50(3,:)
	outmat(1:9,60)=wealthck75(3,:)	
	outmat(1:9,61)=wealthck90(3,:)	
	outmat(1:9,62)=wealthck10(4,:)
	outmat(1:9,63)=wealthck25(4,:)
	outmat(1:9,64)=wealthck50(4,:)
	outmat(1:9,65)=wealthck75(4,:)	
	outmat(1:9,66)=wealthck90(4,:)	
	outmat(1:9,67)=wealthck10(5,:)
	outmat(1:9,68)=wealthck25(5,:)
	outmat(1:9,69)=wealthck50(5,:)
	outmat(1:9,70)=wealthck75(5,:)	
	outmat(1:9,71)=wealthck90(5,:)	
	outmat(1:5,72)=avFTcohort(:)
    outmat(1:6,73)=avinsPO(:)        !mean insurance for low income parents
	
    !Write output matrix to CSV file
    open(unit=outcsv,file='output/outresults.csv',action="write",status="replace")
    write(outcsv,"(73(A,',',:))") 'age','medconp', 'medconk','avincp','avinck', &
    							  'avFT','avFT_sick','avFT_well','avPT','avPT_sick','avPT_well',&
								  'avinf','avins','avmar','avmu','avSIp','avSIk',&
                                  'avhlthy','kapp25','kapp50','kapp75',&
								  'wealthc1p10','wealthc1p25','wealthc1p50','wealthc1p75','wealthc1p90',&
								  'wealthc2p10','wealthc2p25','wealthc2p50','wealthc2p75','wealthc2p90',&
								  'wealthc3p10','wealthc3p25','wealthc3p50','wealthc3p75','wealthc3p90',&
								  'wealthc4p10','wealthc4p25','wealthc4p50','wealthc4p75','wealthc4p90',&
								  'wealthc5p10','wealthc5p25','wealthc5p50','wealthc5p75','wealthc5p90',&
								  'wealthc1k10','wealthc1k25','wealthc1k50','wealthc1k75','wealthc1k90',&
								  'wealthc2k10','wealthc2k25','wealthc2k50','wealthc2k75','wealthc2k90',&
								  'wealthc3k10','wealthc3k25','wealthc3k50','wealthc3k75','wealthc3k90',&
								  'wealthc4k10','wealthc4k25','wealthc4k50','wealthc4k75','wealthc4k90',&
								  'wealthc5k10','wealthc5k25','wealthc5k50','wealthc5k75','wealthc5k90',&
								  'avFTcohort','avinsPO'						  
    do i=1,nt+1
        write(outcsv,"(74(f0.12,',',:))") outmat(i,:)
    end do
    close(outcsv)
endif

	
end subroutine writesim
!----------------------

end module simwrite
